#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>

#define BIN_SIZE 100

#define MAXFLDS 200     /* maximum possible number of fields */
#define MAXFLDSIZE 32   /* longest possible field + 1 = 31 byte field */

#define RMSE 1
#define RSQ 2


void usage(char *name){
	printf("\nUsage:\n");
	printf("%s datafile dvStartCol dvEndCol [optional arguments]\n\nOptional arguments:\n",name);
	printf("  -skip int\tnumber of lines atop the datafile to skip (default is 1)\n");
	printf("  -rmse [float]\tbest fit of predictions to data, expressed as a minimum RMSE value\n");
	printf("  -rsq [float]\tbest fit of predictions to data, expressed as a maximum R^2 value\n");
	printf("  -dvMin float\tminimum value of hypothetical values for all columns (default is 0)\n");
	printf("  -dvMax float\tmaximum value of hypothetical values for all columns (default is 1)\n");
	printf("  -c int\tconvergence criterion x100 (default is 100)\n");
	//printf("  -d delimiters\tdelimiter characters separating columns (default is a tab)\n");
}

void parse(char *record, char *delim, char arr[][MAXFLDSIZE],int *fldcnt){
    char*p=strtok(record,delim);
    int fld=0;
    while(p){
        strcpy(arr[fld],p);
		fld++;
		p=strtok('\0',delim);
	}
	*fldcnt=fld;
}

float frand(float min,float range){return min+range*(float)rand()/(float)RAND_MAX;}

float rsq(int n,float *y,float sumy,float dotyy,float *x,float sumx,float dotxx){
	int i;
	float dotxy=0.0f;
	for(i=0;i<n;i++)dotxy=dotxy+x[i]*y[i];
	//printf("%f %f %f %f %f\n",sumy,dotyy,sumx,dotxx,dotxy);
	return pow(n*dotxy-sumx*sumy, 2) / ((n*dotxx-pow(sumx,2)) * (n*dotyy-pow(sumy,2)));
}

int main(int argc, char **argv) {
	char tmp[1024]={0x0};
	int fldcnt=0;
	char arr[MAXFLDS][MAXFLDSIZE]={0x0};
	int i,j,k;
	FILE *in=fopen(argv[1],"r");
	if(in==NULL){
		printf("Unable to open the file.\nPlease specify a valid path for a file with model predictions.\n");
		exit(EXIT_FAILURE);
	}
	if(argc<4){
		usage(argv[0]);
		fgets(tmp,sizeof(tmp),in);
		parse(tmp,"\t, ",arr,&fldcnt);
		printf("\nFirst line of <%s>:\n", argv[1]);
		for(i=0;i<fldcnt;i++)printf("\tColumn %2d: %s\n",i,arr[i]);
		printf("\n");
	}else{
		int startField=atoi(argv[2]);
		int endField=atoi(argv[3]);
		int fields=endField-startField+1;
		float dataset[fields];
		float dvMin=0.0f,dvMax=1.0f,bestRMSE=-1,bestRSQ=-1,bestSSE=FLT_MAX;
		int stableBinsCriterion=100,skipLines=1;
		int gofMeasures=0;
		char *err;
		double value;
		for(i=4;i<argc;i++){
			err="err";
			if(i+1<argc)value=strtod(argv[i+1],&err);
			if(strcmp(argv[i],"-dvMin")==0){
				if(*err==0){dvMin=(float)value;}
				else{printf("No value for parameter -dvMin\n");exit(0);}
			}else if(strcmp(argv[i],"-dvMax")==0){
				if(*err==0){dvMax=(float)value;}
				else{printf("No value for parameter -dvMin\n");exit(0);}
			}else if(strcmp(argv[i],"-c")==0){
				if(*err==0){stableBinsCriterion=(int)value;}
				else{printf("No value for parameter -s\n");exit(0);}
			}else if(strcmp(argv[i],"-skip")==0){
				if(*err==0){skipLines=(int)value;}
				else{printf("No value for parameter -skip\n");exit(0);}
			}else if(strcmp(argv[i],"-rmse")==0){
				gofMeasures=gofMeasures|RMSE;
				if(*err==0){bestRMSE=value;bestSSE=bestRMSE*bestRMSE*fields;}
			}else if(strcmp(argv[i],"-rsq")==0){
				gofMeasures=gofMeasures|RSQ;
				if(*err==0){bestRSQ=value;}
			}
		}
		if(!gofMeasures)gofMeasures=RMSE;
		//get num of lines
		int recordcnt=0;while(fgets(tmp,sizeof(tmp),in)!=0)recordcnt++;
		fclose(in);
		recordcnt=recordcnt-skipLines;
		printf("total number of records: %d\n", recordcnt);
		printf("number of DV columns: %d\n", endField-startField+1);
		//create predictions matrix
		float **recs=(float**) malloc(recordcnt*sizeof(float*));//[recordcnt][fields];
		float sse,newRSQ,summ[recordcnt],dotmm[recordcnt],sumx=0,dotxx=0;
		in=fopen(argv[1],"r");
		printf("reading file...\n");
		for(j=0;j<skipLines;j++)fgets(tmp,sizeof(tmp),in); //get rid of header
		for(j=0;fgets(tmp,sizeof(tmp),in)!=0;j++){
			recs[j]=(float*) malloc(fields*sizeof(float));
			parse(tmp,"\t, ",arr,&fldcnt);
			for(i=startField;i<=endField;i++)
				recs[j][i-startField]=atof(arr[i]);
		}
		//cache summ and dotmm
		if(gofMeasures&RSQ){
			for(j=0;j<recordcnt;j++){ 
				summ[j]=0.0f;
				dotmm[j]=0.0f;
				for(i=0;i<fields;i++){
					summ[j]=summ[j]+recs[j][i];
					dotmm[j]=dotmm[j]+recs[j][i]*recs[j][i];
				}
			}
		}
		printf("finished reading file.\n");
		//if no gof values are specified, calculate based on observed data
		if((gofMeasures&RMSE && bestRMSE==-1) || (gofMeasures&RSQ && bestRSQ==-1)){
			printf("Enter observed data values (space-separated): ");
			fgets(tmp,1024,stdin);
			parse(tmp,"\t, ",arr,&fldcnt);
			if(fields!=fldcnt){printf("Incorrect number of datapoints.\n");exit(0);}
			for(i=0;i<fields;i++)dataset[i]=atof(arr[i]);
			if(gofMeasures&RMSE){
				for(j=0;j<recordcnt;j++){
					sse=0;
					for(i=0;i<fields;i++)sse=sse+pow(dataset[i]-recs[j][i],2);
					if(sse<bestSSE)bestSSE=sse;
				}
				bestRMSE=sqrt(bestSSE/fields);
				printf("Best fit: RMSE = %f\n",bestRMSE);
			}
			if(gofMeasures&RSQ){
				for(i=0;i<fields;i++){
					sumx=sumx+dataset[i];
					dotxx=dotxx+dataset[i]*dataset[i];
				}
				for(j=0;j<recordcnt;j++){
					newRSQ=rsq(fields,recs[j],summ[j],dotmm[j],dataset,sumx,dotxx);
					if(newRSQ>bestRSQ)bestRSQ=newRSQ;
				}
				printf("Best fit: r^2 = %f\n",bestRSQ);
			}
		}
		//mfa
		srand(time(NULL));
		int B=0,stableBins=0;
		float p=0,lastp=0,dvRange=dvMax-dvMin;
		char pstr[10];
		for(k=1;stableBins<stableBinsCriterion;k++){
			//make hypothetical dataset
			for(i=0;i<fields;i++)dataset[i]=frand(dvMin,dvRange);
			if(gofMeasures&RSQ){
				sumx=0,dotxx=0;
				for(i=0;i<fields;i++){
					sumx=sumx+dataset[i];
					dotxx=dotxx+dataset[i]*dataset[i];
				}
			}
			//check for best fit to hypothetical dataset
			for(j=0;j<recordcnt;j++){
				if(gofMeasures&RMSE){sse=0;for(i=0;i<fields;i++)sse=sse+pow(dataset[i]-recs[j][i],2);}
				if(gofMeasures&RSQ){newRSQ=rsq(fields,recs[j],summ[j],dotmm[j],dataset,sumx,dotxx);}
				if((!gofMeasures&RMSE || sse<=bestSSE) && (!gofMeasures&RSQ || newRSQ>=bestRSQ)){B++;break;}
			}
			//check stability and report MFA
			if(k%BIN_SIZE==0){
				if((float)B/k<.01) sprintf(pstr,"%.03f",(float)B/k);
				else sprintf(pstr,"%.02f",(float)B/k);
				p=atof(pstr);
				if(p==lastp && k>BIN_SIZE)
					stableBins++;
				else {
					stableBins=0;
					lastp=p;
				}
				printf("%d \tφ(",k);
				if(gofMeasures&RMSE)printf("rmse=%.03f,",bestRMSE);
				if(gofMeasures&RSQ)printf("rsq=%.03f,",bestRSQ);
				printf("%d)",stableBins*BIN_SIZE);
				if(p)printf("=%s\n",pstr);
				else printf("<0.001\n");
			}
		}
	}
	fclose(in);
    return 0;	
}
